Comparación de distintos modelos

para la predicción de Series Temporales

Un enfoque híbrido en Python y R

0 - Nota aclaratoria

Toda la narrativa y explicaciones encontradas a lo largo del trabajo están hechas en español para que su correción sea más sencilla. Pero por consistencia con mis hábitos de programación, todo el código hecho en Python o R está en inglés (variables, funciones, comentarios dentro del código).

Mi intención es la de traducir todo el contenido en castellano al inglés y así poder subir el código y el resto de las explicaciones como un kernel a Kaggle como se recomendaba en la guía del TFM. Cuando tenga todo el contenido traducido al inglés lo subiré a Kaggle con el siguiente perfil: https://www.kaggle.com/avallvaq/kernels con el título - Comparison of time series forecast models. A hybrid approach in Python and R - En este apartado incluiré la nota: This analysis is part of the final project for the Master's Degree: "Master de BigData – UNED"

1 - Introducción

Se trada de un subconjunto de datos numéricos, entre los que existe una relación temporal.

Las series temporales necesitan su propia metodología de análisis, ya que la componente temporal añade un orden de complejidad y autocorrelación entre los datos de una misma variable.

Igual que en datos estacionarios se buscan patrones "ocultos" en los datos, ya sea entre variables del data set o variables creadas mediante proceso de data engineering. Para las series temporales se buscan la existencia de estacionalidad, análisis de tendencias o fluctuaciones periódicas (ciclos).

Existen diversos algoritmos que se pueden utilizar para analizar series temporales. Desde regresiones lineales sencillas, para hacer análisis de tendencia más sencillos, regresión mediante árboles de decisión, Random Forest o Gradient Boosting Methods o modelos ARIMA (AutoRegressive Integrated Moving Average). También existen modelos más complejos que pueden llegar a capturar las relaciones no lineales de series temporales complejas, como pueden ser las redes neuronales LSTM (Long-short Term Memory Networks), TDNN (Time Delay Neural Networks).

1.1 - Dataset de Kaggle: Hourly Energy Consumption

El dataset puede descargarse de Kaggle en la siguiente URL: https://www.kaggle.com/robikscube/hourly-energy-consumption#NI_hourly.csv

Se trata de una serie temporal univariante, es decir, que sólo aporta datos históricos de una sóla variable. En este caso se trata de una estimación del consumo eléctrico en Mega Watios (MW) donado por la compañia PJM Interconection LLC. Se trata de una compañía de distribución eléctrica regional en los Estados Unidos. Se encarga de distribuir electricidad entre varios estados de la zona noreste y medio este del país (algunos de ellos, Delaware, Illinois, Indiana, Kentucky, Maryland, Michigan, New Jersey, Ohio, Pennsylvania, Tennessee, Virginia, etc.)

1.2 - Análisis

El trabajo se ha realizado utilizando tanto Python como R. En Python se ha llevado a cabo la mayoría del EDA aprovechando la flexibilidad y potencia de la librería pandas, la mayoríad de los gráficos se renderizan en plotly, ya que la interactividad ha resultado muy útil en el contexto del análisis de series. También se ha recurrido a matplotlib y seaborn para algunas de las visualizaciones. El análisis de la serie temporal se ha llevado a cabo con la librería statsmodel y la librería Prophet de Facebook.

R es conocido por su capacidad a la hora de hacer análisis y predicción de series temporales con los objetos ts. Para el análisis de los datos se ha optado por utilizar la librería forecast, que a su vez se apoya en otros paquetes vistos durante el curso como ggplot.

El trabajo se ha desarrollado de la siguiente forma:

  • Sección 2. Ingestión de datos y limpieza (eliminación de valores duplicados, creación del índice temporal, etc)

  • Sección 2.1. EDA: Búsqueda de patrones estacionales, descomposición de la serie en tendecia, estacionalidad y residuos

  • Sección 3. Comprobación formal de la estacionalidad de la serie. Contraste de la hipótesis de raíz unitaria (método KPPS)

  • Sección 4. División simple en conjuntos de entrenamiento y validación

  • Secciones 5, 6, 7, 8 y 9. Predicciones

  • Sección 10. Tabla con los valores MAPE de las predicciones

1.3 - Métodos de predicción utilizados

En Python:

  • I - Holt-Winters con estacionalidad

  • II - Prophet

En R:

  • III - Métodos sencillos en R (media, naïve, etc)
  • IV - Regresión armónica dinámica con series de Fourier
  • V - TBATS
In [1]:
import numpy as np
import pandas as pd
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
# import pandas_profiling
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_squared_error, mean_absolute_error
import statsmodels.api as sm
import xgboost as xgb
from fbprophet import Prophet
from fbprophet.diagnostics import cross_validation

# graphic visualization packages
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import seaborn as sns
import chart_studio.plotly as py
import plotly.graph_objects as go
import plotly.io as pio
import plotly.express as px
import plotly.figure_factory as ff

%matplotlib inline

1.5 - Paquetes necesarios en R

# Packages needed for the R environment:
library(tidyverse)
library(ggfortify)
library(lubridate)
library(forecast)
library(tseries)
library(urca)
In [2]:
# Read the data into a pandas DataFrame
pjme = pd.read_csv("data/PJME_hourly.csv", parse_dates=[0], encoding='UTF-8')
pjme.Datetime = pd.to_datetime(pjme.Datetime)
pjme.rename(columns={'PJME_MW':'demand_in_MW'}, inplace=True)

# Drop possible duplicate values and order the data by date and time
pjme.drop_duplicates(subset='Datetime', keep='last', inplace=True)
pjme.sort_values(by=['Datetime'], axis=0, ascending=True, inplace=True)
# Reindex after the dataframe has been sorted
pjme.reset_index(drop=True, inplace=True)


# Now, let's check if the dataframe have some missing measurments
pjme = pjme.set_index('Datetime')
print(f'Datetime freq is set to: {pjme.index.freq}, i.e. the dataset is not continous') # Datatime index frequency is set to None
                                                     # i.e. the time series is not cointinous

# We can measure the number of missing measurements
date_range = pd.date_range(start=min(pjme.index), 
                           end=max(pjme.index), 
                           freq='H')
print(f'The number of missing measurements in the set is: {(len(date_range)-len(pjme))}:')

# This will append the previously missing datetimes, and create null values in our target variable
pjme = pjme.reindex(date_range)
# Fill in the blanks with values that lie on a linear curve between existing data points
pjme['demand_in_MW'].interpolate(method='linear', inplace=True)
# Now the dataset is continously populated
print(f'The pjme.index.freq is now: {pjme.index.freq}, indicating that we no longer have missing instances')
# Lets copy the Datetime index again into a columns, we are going to need for further calculations
pjme['Datetime'] = pjme.index
pjme.Datetime = pd.to_datetime(pjme.Datetime)

print("\n")

# Display the first and last rows
display(pjme.head())
display(pjme.tail())

print(f"Length of the dataset: {len(pjme):,}\n\n\
Number of N/A values:\n{pd.isna(pjme).sum()}\n\
Number of missing values:\n{pd.isnull(pjme).sum()}")
Datetime freq is set to: None, i.e. the dataset is not continous
The number of missing measurements in the set is: 30:
The pjme.index.freq is now: <Hour>, indicating that we no longer have missing instances


demand_in_MW Datetime
2002-01-01 01:00:00 30393.0 2002-01-01 01:00:00
2002-01-01 02:00:00 29265.0 2002-01-01 02:00:00
2002-01-01 03:00:00 28357.0 2002-01-01 03:00:00
2002-01-01 04:00:00 27899.0 2002-01-01 04:00:00
2002-01-01 05:00:00 28057.0 2002-01-01 05:00:00
demand_in_MW Datetime
2018-08-02 20:00:00 44057.0 2018-08-02 20:00:00
2018-08-02 21:00:00 43256.0 2018-08-02 21:00:00
2018-08-02 22:00:00 41552.0 2018-08-02 22:00:00
2018-08-02 23:00:00 38500.0 2018-08-02 23:00:00
2018-08-03 00:00:00 35486.0 2018-08-03 00:00:00
Length of the dataset: 145,392

Number of N/A values:
demand_in_MW    0
Datetime        0
dtype: int64
Number of missing values:
demand_in_MW    0
Datetime        0
dtype: int64

2.1 - Análisis exploratorio de los datos (EDA)

In [3]:
# Let's create some seasonal features
def create_features(df, label=None):
    """
    Function to create several time-range features, day of the week, day of the year, day of month, 
    """
    
    df['date'] = df.Datetime
    df['hour'] = df['date'].dt.hour
    df['dayofweek'] = df['date'].dt.dayofweek
    df['quarter'] = df['date'].dt.quarter
    df['month'] = df['date'].dt.month
    df['year'] = df['date'].dt.year
    df['dayofyear'] = df['date'].dt.dayofyear
    df['dayofmonth'] = df['date'].dt.day
    df['weekofyear'] = df['date'].dt.weekofyear   
    # Create the season feature
    conditions = [
            (df['dayofyear'] < 79),
            (df['dayofyear'] >= 356),
            (df['dayofyear'] >= 79)  & (df['dayofyear'] < 172),
            (df['dayofyear'] >= 172) & (df['dayofyear'] < 266),
            (df['dayofyear'] >= 266  & (df['dayofyear'] < 356))]
    choices = [0, 0, 1,2,3]
    df['season'] = np.select(conditions, choices, default=0)
        
    X = df[['date','hour','dayofweek','month', 'quarter', 'season',
           'dayofyear', 'dayofmonth', 'weekofyear', 'year']]
    if label:
        y = df[label]
        return pd.concat([X, y], axis=1)
    return X
In [4]:
def map_months_days_to_string(df):
    """
    Simple function to rename month and days of the week, from numeric values to chars.
    """

    df['dayofweek_string'] = df['dayofweek'].map({0:'Monday', 1:'Tuesday', 2:'Wednesday', 3:'Thursday',
                                                  4: 'Friday', 5:'Saturday', 6:'Sunday'})
    df['season_string']    = df['season'].map({0:'winter', 1:'spring', 2:'summer', 3:'fall'})
    
    return df[['date','hour','dayofweek','dayofweek_string', 'month', 'quarter', 'season',
               'season_string', 'dayofyear', 'dayofmonth', 'weekofyear', 'year', 'demand_in_MW']]
In [5]:
# I use the features_target dataframes for the EDA
features_target = create_features(pjme, 'demand_in_MW')
features_target = map_months_days_to_string(features_target)

# I just mapped numerical categories to strings of the actual days of the week and seasons
# for a better interpretability of the group by aggregation computations in the visual EDA section
display(features_target.head())
display(features_target.tail())
date hour dayofweek dayofweek_string month quarter season season_string dayofyear dayofmonth weekofyear year demand_in_MW
2002-01-01 01:00:00 2002-01-01 01:00:00 1 1 Tuesday 1 1 0 winter 1 1 1 2002 30393.0
2002-01-01 02:00:00 2002-01-01 02:00:00 2 1 Tuesday 1 1 0 winter 1 1 1 2002 29265.0
2002-01-01 03:00:00 2002-01-01 03:00:00 3 1 Tuesday 1 1 0 winter 1 1 1 2002 28357.0
2002-01-01 04:00:00 2002-01-01 04:00:00 4 1 Tuesday 1 1 0 winter 1 1 1 2002 27899.0
2002-01-01 05:00:00 2002-01-01 05:00:00 5 1 Tuesday 1 1 0 winter 1 1 1 2002 28057.0
date hour dayofweek dayofweek_string month quarter season season_string dayofyear dayofmonth weekofyear year demand_in_MW
2018-08-02 20:00:00 2018-08-02 20:00:00 20 3 Thursday 8 3 2 summer 214 2 31 2018 44057.0
2018-08-02 21:00:00 2018-08-02 21:00:00 21 3 Thursday 8 3 2 summer 214 2 31 2018 43256.0
2018-08-02 22:00:00 2018-08-02 22:00:00 22 3 Thursday 8 3 2 summer 214 2 31 2018 41552.0
2018-08-02 23:00:00 2018-08-02 23:00:00 23 3 Thursday 8 3 2 summer 214 2 31 2018 38500.0
2018-08-03 00:00:00 2018-08-03 00:00:00 0 4 Friday 8 3 2 summer 215 3 31 2018 35486.0
In [6]:
fig = go.Figure(layout=go.Layout(title=go.layout.Title(
    text=f'Power Demand (MW) over time [{min(pjme.year)} - {max(pjme.year)}]')))
fig.add_trace(go.Scatter(x=pjme.date, y=pjme.demand_in_MW, name='Demand in MW', line_color='chocolate'))
fig.update_traces(line=dict(width=0.5))
fig.update_xaxes(title_text='Time')
fig.update_yaxes(title_text='Power Consumption (MW)')
fig.update_layout(showlegend=True, width=1000)

fig.show()
In [7]:
# Compute the correlation matrix
corr = features_target[['demand_in_MW', 'hour', 'dayofweek', 'month', 'season']].corr()

# Generate a mask for the upper triangle
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(8, 6))

# Generate a custom diverging colormap
colormap = sns.diverging_palette(220, 10, as_cmap=True)

# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr, mask=mask, cmap=colormap, center=0, vmin=-1.0, vmax=1.0,
            square=True, linewidths=.5, cbar_kws={"shrink": 1}, annot=True)
Out[7]:
<matplotlib.axes._subplots.AxesSubplot at 0x1fe83832048>

Podemos ver, como la hora del día tiene una correlación positiva. Por lo que en principio el consumo subirá en horas más entrado el día. A su vez también podemos ver una ligera correlación negativa del consumo eléctrico con el día de la semana. Teniendo en cuenta que hemos codificado los fines de semana como los días 5 y 6 de la semana, nos lleva a esperar una menor demanda durante el fin de semana.

2.1.1 - En busca de patrones temporales

In [8]:
# Let's aggregate the data by hour of the day during the week
features_target_aggregated = features_target.groupby(['hour', 'dayofweek_string'], as_index=False).agg({'demand_in_MW':'median'})

#plot
fig = px.line(features_target_aggregated, x='hour', y='demand_in_MW', 
              color='dayofweek_string',
              title='Median Hourly Power Demand per Weekday')
fig.update_traces(mode='lines+markers')
fig.update_layout(xaxis_title='Hour',
                  yaxis_title='Energy Demand [MW]')
fig.for_each_trace(
    lambda trace: trace.update(name=trace.name.replace("dayofweek_string=", "")),
)
fig.show()

Podemos observar como los fines de semana la demanda de electricidad es alrededor de un 10% menor que entre semana. A priori sería entendible plantear la hipótesis del consumo eléctrico durante la semana con el consumo que generan los puestos de trabajo.

In [9]:
# Let's aggregate by hour of the day during different seasons
features_target_aggregated_by_season = features_target.groupby(['hour', 'season_string'], as_index=False)\
                                        .agg({'demand_in_MW':'median'})

# plot
color_map_season = {'winter':'mediumpurple', 'spring': 'forestgreen', 'summer': 'gold', 'fall':'chocolate'}
fig = px.line(features_target_aggregated_by_season, x='hour', y='demand_in_MW', 
              color='season_string', title='Median Hourly Power Demand per Season', color_discrete_map=color_map_season)
fig.update_traces(mode='lines+markers')
fig.update_layout(xaxis_title='Hour',
                  yaxis_title='Energy Demand [MW]')
fig.for_each_trace(
    lambda trace: trace.update(name=trace.name.replace("season_string=", "")),
)
fig.show()

Aunque el coeficiente de correlación de Pearson era cercano a cero. Era esperado encontrar una estacionalidad alta del consumo eléctrico. Vemos como la demanda está muy estrechamente ligado al verano y en menor medida al invierno.

2.1.2 - Descomposición de la Serie Temporal

Uno de los aspectos más interesantes de las series temporales, es el poder descomponer la serie en diferentes componentes. Una componente de soporte principal de la tendencia, una componente periódica, que es la que caracteriza los ciclos estacionales.

Como en este ejemplo la componente estacional es constante, como hemos podido ver en la Figura 1.

In [10]:
from statsmodels.tsa.seasonal import seasonal_decompose
# NOTE, that seasonal_decompose needs to have the dataframe with a Datetime object as index
series = features_target[['demand_in_MW']]
frequency = 24*365

# decomposing the time-series, with the frequency being 24 hours per 365 days
decomposed = seasonal_decompose(series, model='additive', freq=frequency)
In [11]:
def plot_decompositions(decompositions, titles, line_widths):
    " Plotting the different elements constituting the time-series"
    
    for data, title, lw in zip(decompositions, titles, line_widths):
        # draw a line plot of the data with plotly express
        fig = px.line(data, y='demand_in_MW',
                      title=title, height=300)
        # adjust line width
        fig.update_traces(line=dict(width=lw), line_color='chocolate')
        
        fig.update_layout(xaxis=dict(showticklabels=False, linewidth=1),
                          yaxis=dict(title=''), 
                          margin=go.layout.Margin(l=40, r=40, b=0, t=40, pad=0),)
        if title=='Trend':
            fig.update_yaxes(range=[10000,65000])
        
        fig.update_layout(xaxis_title="Years")
        fig.show()

# calling the function 
plot_decompositions(decompositions=[decomposed.trend, decomposed.seasonal, decomposed.resid],
                    titles=['Trend', 'Seasonality', 'Residuals'],
                    line_widths=[2, 0.05, 0.5])

La razón del modelo aditivo, es porque se trata de un modelo lineal, en el que la componente de tendencia (trend) tiende a ser constante, la componente estacional (seasonal) tiene una periodicidad (amplitud y frecuencia constantes).

3 - Comprobación de la estacionariedad de la serie (hipótesis de raíz unitaria) - R

A priori, debido a que la tendencia de la serie temporal es visualmente muy "plana", podríamos pensar que se trata de una serie estacionaria. Pero aún así, vamos a realizar una comprobación más formal de ello utilizando el método KPSS para comprobar la hipótesis de la raíz unitaria.

Hay varios métodos a disposicón para comprobar la estacionalidad de una serie temporal: -Métodos cualitativos visuales (gráfica de los datos, gráficas de autocorrelación, estadísticos móviles, etc.) -Métodos cuantitativos como el Alternative Dickey-Fuller Test (ADF) o el método KPSS

Primer punto donde introducimos código en R. Vamos a proceder a evaluar la estacionariedad de los datos, utilizando el método KPSS incluido en el paquete forecast() de R.

# First, we import the data into the R workspace
pjme = read.csv("data/PJME_hourly.csv")
pjme[1] <- lapply(pjme[1], as.character)
pjme[,1]  <-  parse_date_time(pjme[,1], "ymd HMS")
pjme <- distinct(pjme,Datetime, .keep_all = TRUE)
pjme <- arrange(pjme,Datetime)
pjme <- pjme[order(pjme$Datetime),]

# Train/Test split. Train: From 1st Jan 2016 to 30th Nov 2017. Test: From 1st Dec 2017 onwards
train <- pjme %>% 
  filter(Datetime < ymd_hms("2017-12-01 00:00:00")) %>% 
  filter(Datetime > ymd_hms("2015-12-31 23:00:00"))
test  <- pjme %>% 
  filter(Datetime > ymd_hms("2017-11-30 23:00:00"))

# Create the msts() object. Which is basically a Time Series object that accepts multiple seasonality
train_ts     <- msts(train$PJME_MW, start = c(2016,0), ts.frequency = 24*365.25, seasonal.periods = c(24,24*365.25))
# Since we analyzed the data and reach to the conclusion that weekly seasonality is not as strong as daily and yearly,
# I did not included it in the seasonal.periods vector, thus c(24, 24*365.25).



train_ts %>% diff() %>% ur.kpss() %>% summary()

ndiffs(train_ts) # We should take a first difference for the data

train_ts %>% diff(lag=24) %>% ndiffs() # No seasonal difference.

Pasamos a realizar las pruebas de raíz unitaria con la función KPSS del paquete urca.

# ====================
# Stationarity check
# ====================
>>> train_ts %>% ur.kpss() %>% summary() # Non-differenced data
####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 14 lags. 

Value of test-statistic is: 2.2087 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

2.21 es un valor demasiado alto del estadístico, y por lo tanto no podemos descartar la hipótesis nula de no haber una raíz unitaria en la serie temporal. Procedo a diferenciar la serie una vez más:

>>> train_ts %>% diff() %>% ur.kpss() %>% summary()
####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 14 lags. 

Value of test-statistic is: 8e-04 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

>>> ndiffs(train_ts) # We should take ONE first difference for the data
1

Tras realizar una diferenciación, los datos son estacionarios. Este es un punto muy importante del análisis, ya que el parámetro d en cualquier modelo de la familia ARIMA se decide bajo esta premisa, el número de veces que hay que diferenciar la serie hasta estacionarizarla es el valor de d.

Hacemos lo mimso con la diferencia estaciona y vemos que el resultado es nulo, por lo que no sería necesario diferenciar para hacer estacionaria la serie temporal. D=0. Es importante también tener en cuenta que el valor de D se puede ver a simple vista ya que la estacionalidad diaria (lag=24) está muy marcada y se ve que es fácilmente predecible.

>>> train_ts %>% diff(lag=24) %>% ndiffs() # No daily seasonal difference needed
0
>>> train_ts %>% diff(lag=168) %>% ndiffs() # No weekly seasonal difference needed
0
In [12]:
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Let's plot the hourle autocorrelation.
fig, ax = plt.subplots(2, figsize=(12,6))
ax[0] = plot_acf(features_target.demand_in_MW.diff().dropna(), lags=168, ax=ax[0])
ax[1] = plot_pacf(features_target.demand_in_MW.diff().dropna(), lags=168, ax=ax[1])

Vemos en la gráfica de autocorrelación como la estacionalidad diaria es muy importante, ya que cada 24 horas tenemos un pico en la gráfica de autocorrelación.

Vamos a utilizar ahora la función mstl() del paquete forecast de R. Esta función nos permite descomponer la serie en distintas estacionalidades. Esto va a permitir ver además de la componente diaria, también las componentes de variación semanales y anuales.

# ============================================
# STL decomposition for multiple seasonality
# ============================================
pjme_ts %>% mstl() %>% 
  autoplot() + xlab("Year")

La gráfica se lee de la siguiente forma:

  • Señal en crudo
  • Tendencia
  • Componente estacional diaria (frecuencia m=24)
  • Componente estacional semanal (frecuencia m=24x7=168)
  • Componente estacional anual (frecuencia m=24x365.25=8766)

Es importante fijarse en la escala vertical de las gráficas. Especialmente para ver como la estacionalidad semanal es del orden un 50% más débil que la estacionalidad diaria y anual.

Para el proyecto, vamos a reducir los datos a un periodo de dos años y un trimestre. Debido a la frecuencia horaria de este data-set y de la estacionalidad diaria y anual, un periodo de dos años puede contener los suficientes datos como para que nuestros modelos realicen sus estimaciones, además, reducimos el coste computacional de los modelos de predicción.

En principio y por temas de tiempo de ejecución de los distintos modelos. He obtado por una evaluación simple, con una división de los datos en entrenamiento y test del 66% 34% de la longitud de la serie. Ver Anexo para ver más información acerca de la validación cruzada de series temporales.

In [13]:
def ts_train_test_split(df, init_date, end_date, cutoff_date):
    """
    Function to return a simple train/test split of the given a univariate time-series, within dates provided by the user.
    The date should be given as a string in ISO format YYYY-MM-DD.
    """
    
    start = init_date
    end   = end_date
    print(f'The first date time point is: {start}')
    print(f'The last date time point is: {max(df.index)}')
    
    range_df = len(pd.date_range(start=start,end=end, freq='Y'))
    print(f"The difference is: {range_df} years")
    print(f'The training range has to be {range_df*0.66} years')
        
    CUTOFF_DATE = pd.to_datetime(cutoff_date)
    train = df.loc[(df.index < CUTOFF_DATE) & (df.index > init_date)].copy(deep=True)
    test = df.loc[df.index >= CUTOFF_DATE].copy(deep=True)

    return train, test

train, test = ts_train_test_split(features_target, '2016-01-01 00:00:00', max(features_target.index), '2017-12-01')
display(train.tail())
display(test.head())

y_train = train['demand_in_MW']
y_test  = test['demand_in_MW']
display(y_test)
The first date time point is: 2016-01-01 00:00:00
The last date time point is: 2018-08-03 00:00:00
The difference is: 2 years
The training range has to be 1.32 years
date hour dayofweek dayofweek_string month quarter season season_string dayofyear dayofmonth weekofyear year demand_in_MW
2017-11-30 19:00:00 2017-11-30 19:00:00 19 3 Thursday 11 4 3 fall 334 30 48 2017 33825.0
2017-11-30 20:00:00 2017-11-30 20:00:00 20 3 Thursday 11 4 3 fall 334 30 48 2017 33377.0
2017-11-30 21:00:00 2017-11-30 21:00:00 21 3 Thursday 11 4 3 fall 334 30 48 2017 32491.0
2017-11-30 22:00:00 2017-11-30 22:00:00 22 3 Thursday 11 4 3 fall 334 30 48 2017 30845.0
2017-11-30 23:00:00 2017-11-30 23:00:00 23 3 Thursday 11 4 3 fall 334 30 48 2017 28581.0
date hour dayofweek dayofweek_string month quarter season season_string dayofyear dayofmonth weekofyear year demand_in_MW
2017-12-01 00:00:00 2017-12-01 00:00:00 0 4 Friday 12 4 3 fall 335 1 48 2017 26364.0
2017-12-01 01:00:00 2017-12-01 01:00:00 1 4 Friday 12 4 3 fall 335 1 48 2017 24804.0
2017-12-01 02:00:00 2017-12-01 02:00:00 2 4 Friday 12 4 3 fall 335 1 48 2017 24003.0
2017-12-01 03:00:00 2017-12-01 03:00:00 3 4 Friday 12 4 3 fall 335 1 48 2017 23780.0
2017-12-01 04:00:00 2017-12-01 04:00:00 4 4 Friday 12 4 3 fall 335 1 48 2017 23891.0
2017-12-01 00:00:00    26364.0
2017-12-01 01:00:00    24804.0
2017-12-01 02:00:00    24003.0
2017-12-01 03:00:00    23780.0
2017-12-01 04:00:00    23891.0
                        ...   
2018-08-02 20:00:00    44057.0
2018-08-02 21:00:00    43256.0
2018-08-02 22:00:00    41552.0
2018-08-02 23:00:00    38500.0
2018-08-03 00:00:00    35486.0
Freq: H, Name: demand_in_MW, Length: 5881, dtype: float64
In [14]:
# Plot Train/Test
fig = go.Figure(layout=go.Layout(title=go.layout.Title(
    text=f'Power Demand (MW) over time [{min(train.year)} - {max(test.year)}]')))
fig.add_trace(go.Scatter(x=train.date, y=train.demand_in_MW, name='Training', line_color='chocolate'))
fig.add_trace(go.Scatter(x=test.date, y=test.demand_in_MW, name='Test', line_color='green'))
fig.update_traces(line=dict(width=0.5))
fig.update_xaxes(title_text='Time')
fig.update_yaxes(title_text='Power Consumption (MW)')
fig.update_layout(showlegend=True, width=1000)

fig.show()

5 - Método I: Holt-Winters con triple suavizamiento exponencial

El algoritmo de Holt-Winters es un algoritmo basado en medias móviles muy utilizado para la predicción de series temporales. Debido a su simplicidad e interpretabilidad, vamos a utilizarlo como algoritmo base con el que después comparar los siguientes.

La mayor limitación de este algoritmo es la capacidad de tener en cuenta sólo una estacionalidad. Nos decantamos por tomar el periodo anual, seasonal_periods=24*365

In [15]:
import statsmodels.api as sm
# exponential smoothing only takes into consideration patterns in the target variable
# so we discard the other features

# fit & predict
holt_winter = sm.tsa.ExponentialSmoothing(y_train, seasonal='add', seasonal_periods=24*365).fit()
y_hat_holt_winter = holt_winter.forecast(len(y_test))
C:\Users\avv1g\Anaconda3\lib\site-packages\statsmodels\tsa\holtwinters.py:712: ConvergenceWarning:

Optimization failed to converge. Check mle_retvals.

In [16]:
# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=y_test.index, y=y_test,
                         mode='lines',
                         name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=y_hat_holt_winter.index, y=y_hat_holt_winter,
                         mode='lines', 
                         name='Test - Prediction'))

# adjust layout
fig.update_traces(line=dict(width=0.5))
fig.update_layout(title='Holt-Winter Forecast of Hourly Energy Demand',
                  xaxis_title='Date & Time',
                  yaxis_title='Energy Demand [MW]')

A primera vista el modelo parece obtener una buena aproximación de las series temporales. Pero vamos a evaluarlo de calculando el Mean Absolute Percentage Error (MAPE) entre el conjunto de validación y_test y la predicción y_hat_holt_winter.

In [17]:
def mape(y_true, y_pred):
    """ 
    Mean Absolute Percentage Error
    """
    
    # convert to numpy arrays
    y_true, y_pred = np.array(y_true), np.array(y_pred)
    
    # take the percentage error
    pe = (y_true - y_pred) / y_true
    
    # take the absolute values
    ape = np.abs(pe)
    
    # quantify the performance in a single number
    mape = np.mean(ape)
    
    return f'{mape*100:.2f}%'
In [18]:
mape_hw = mape(y_true=y_test, y_pred=y_hat_holt_winter)
print(f'Our Holt-Winters model has a mean average percentage error (MAPE) of: {mape_hw}')
Our Holt-Winters model has a mean average percentage error (MAPE) of: 13.11%

A continuación se muestran las gráficas de los resultados en una ventana de tiempo menor (semanal).

In [19]:
# interval length
interval = 24 * 7

# intermediary variables for readability
x_true, y_true = y_test.iloc[:interval].index, y_test.iloc[:interval]
x_pred, y_pred = y_hat_holt_winter.iloc[:interval].index, y_hat_holt_winter.iloc[:interval]

# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_true, y=y_true,
                         mode='lines',
                         name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=x_pred, y=y_pred,
                         mode='lines', 
                         name='Test - Prediction'))

# adjust layout
fig.update_traces(line=dict(width=0.9))
fig.update_layout(title=f'Holt-Winter Intra-Day Forecast of First {interval} Hours of Energy Demand',
                  xaxis_title='Date & Time',
                  yaxis_title='Energy Demand [MW]')
fig.show()

# quantify accuracy
print(f'MAPE for interval of the first {interval} hours: {mape(y_true, y_pred)}')
MAPE for interval of the first 168 hours: 5.88%

La misma ventana de tiempo pero para la última semana de predicciones.

In [20]:
# interval length
interval = -24 * 7

# intermediary variables for readability
x_true, y_true = y_test.iloc[interval:].index, y_test.iloc[interval:]
x_pred, y_pred = y_hat_holt_winter.iloc[interval:].index, y_hat_holt_winter.iloc[interval:]

# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_true, y=y_true,
                         mode='lines',
                         name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=x_pred, y=y_pred,
                         mode='lines', 
                         name='Test - Prediction'))

# adjust layout
fig.update_traces(line=dict(width=0.9))
fig.update_layout(title=f'Holt-Winter Intra-Day Forecast of Last {abs(interval)} Hours of Energy Demand',
                  xaxis_title='Date & Time (yyyy/mm/dd hh:MM)',
                  yaxis_title='Energy Demand [MW]')
fig.show()

# quantify accuracy
print(f'MAPE for interval of the last {abs(interval)} hours: {mape(y_true, y_pred)}')
MAPE for interval of the last 168 hours: 14.24%

Holt-Winters MAPE: 13.11%

1st week MAPE: 5.88%

Last week MAPE: 14.24%

5.1 - El problema de la estacionalidad múltiple

El método de Holt-Winters con estacionalidad no está pensado para capturar estacionalidades múltiples. Al definir la frecuencia de muestreo nos vimos obligados a dar solamente un valor (24x365). A continuación vamos a utilizar modelos mas sofisticados que permiten definir estacionalidad múltiple.

6 - Método II: Predicción utilizando la librería Prophet de Facebook

Enlace al white paper de Prophet: https://peerj.com/preprints/3190/

Prophet se define como un método modular de regresión. Al ser modular permite muy fácilmente añadir o quitar parámetros de manera muy intuitiva, ideal para que analistas y data scientist puedan hacer uso de él y sacarle toda su potencia con una curva de aprendizaje muy suave.

La API está disponible tanto en Python como R, y es rápido, completamente automático aunque lo suficientemente flexible para que el usuario final pueda ir ajustando los modelos en función de sus necesidades o de forma iterativa.

La ventaja que tiene Prophet respecto al Método de Holt-Winters implementado en el apartado anterior, es que es capaz de lidiar con estacionalidades múltiples, eventos de un día (Black Friday, Final del Mundial...), fiestas nacionales (como Navidad, Año Nuevo, etc), periodos vacacionales, etc.

In [21]:
from fbprophet import Prophet
from fbprophet.diagnostics import cross_validation

# format data for prophet model using 'ds' and 'y' as per the Python API
train_prophet = y_train.reset_index().rename(columns={'index':'ds', 'demand_in_MW':'y'})

test_prophet = y_test.reset_index().rename(columns={'index':'ds', 'demand_in_MW':'y'})

display(train_prophet.tail())
display(test_prophet.head())
ds y
16794 2017-11-30 19:00:00 33825.0
16795 2017-11-30 20:00:00 33377.0
16796 2017-11-30 21:00:00 32491.0
16797 2017-11-30 22:00:00 30845.0
16798 2017-11-30 23:00:00 28581.0
ds y
0 2017-12-01 00:00:00 26364.0
1 2017-12-01 01:00:00 24804.0
2 2017-12-01 02:00:00 24003.0
3 2017-12-01 03:00:00 23780.0
4 2017-12-01 04:00:00 23891.0
In [22]:
# Conditions for the variable seasonality of the data
def is_spring(ds):
    date = pd.to_datetime(ds)
    return (date.dayofyear >= 79) & (date.dayofyear < 172)

def is_summer(ds):
    date = pd.to_datetime(ds)
    return (date.dayofyear >= 172) & (date.dayofyear < 266)

def is_autumn(ds):
    date = pd.to_datetime(ds)
    return (date.dayofyear >= 266) & (date.dayofyear < 356)

def is_winter(ds):
    date = pd.to_datetime(ds)
    return (date.dayofyear < 79) | (date.dayofyear >= 356)

def is_weekend(ds):
    date = pd.to_datetime(ds)
    return date.day_name in ('Saturday', 'Sunday')

# adding to train set
train_prophet['is_spring'] = train_prophet['ds'].apply(is_spring)
train_prophet['is_summer'] = train_prophet['ds'].apply(is_summer)
train_prophet['is_autumn'] = train_prophet['ds'].apply(is_autumn)
train_prophet['is_winter'] = train_prophet['ds'].apply(is_winter)
train_prophet['is_weekend'] = train_prophet['ds'].apply(is_weekend)
train_prophet['is_weekday'] = ~train_prophet['ds'].apply(is_weekend)

# adding to test set
test_prophet['is_spring'] = test_prophet['ds'].apply(is_spring)
test_prophet['is_summer'] = test_prophet['ds'].apply(is_summer)
test_prophet['is_autumn'] = test_prophet['ds'].apply(is_autumn)
test_prophet['is_winter'] = test_prophet['ds'].apply(is_winter)
test_prophet['is_weekend'] = test_prophet['ds'].apply(is_weekend)
test_prophet['is_weekday'] = ~test_prophet['ds'].apply(is_weekend)

display(train_prophet.tail())
display(test_prophet.head())
ds y is_spring is_summer is_autumn is_winter is_weekend is_weekday
16794 2017-11-30 19:00:00 33825.0 False False True False False True
16795 2017-11-30 20:00:00 33377.0 False False True False False True
16796 2017-11-30 21:00:00 32491.0 False False True False False True
16797 2017-11-30 22:00:00 30845.0 False False True False False True
16798 2017-11-30 23:00:00 28581.0 False False True False False True
ds y is_spring is_summer is_autumn is_winter is_weekend is_weekday
0 2017-12-01 00:00:00 26364.0 False False True False False True
1 2017-12-01 01:00:00 24804.0 False False True False False True
2 2017-12-01 02:00:00 24003.0 False False True False False True
3 2017-12-01 03:00:00 23780.0 False False True False False True
4 2017-12-01 04:00:00 23891.0 False False True False False True
In [23]:
# instantiating the Prophet class with user defined settings
prophet = Prophet(daily_seasonality=False, weekly_seasonality=False,
                  yearly_seasonality=False)

# custom seasonalities to account for conditional variance 
# (more extreme trends in extreme seasons)
prophet.add_seasonality(name='yearly', period=365.25, fourier_order=10)
prophet.add_seasonality(name='weekly_spring', 
                        period=7,
                        fourier_order=5, 
                        condition_name='is_spring')
prophet.add_seasonality(name='weekly_summer', 
                        period=7,
                        fourier_order=5, 
                        condition_name='is_summer')
prophet.add_seasonality(name='weekly_autumn', 
                        period=7,
                        fourier_order=5, 
                        condition_name='is_autumn')
prophet.add_seasonality(name='weekly_winter', 
                        period=7,
                        fourier_order=5, 
                        condition_name='is_winter')
prophet.add_seasonality(name='daily_spring',  
                        period=1,
                        fourier_order=5, 
                        condition_name='is_spring')
prophet.add_seasonality(name='daily_summer',  
                        period=1,
                        fourier_order=5, 
                        condition_name='is_summer')
prophet.add_seasonality(name='daily_autumn',  
                        period=1,
                        fourier_order=5, 
                        condition_name='is_autumn')
prophet.add_seasonality(name='daily_winter',  
                        period=1,
                        fourier_order=5, 
                        condition_name='is_winter')
prophet.add_seasonality(name='daily_weekend',  
                        period=1,
                        fourier_order=5, 
                        condition_name='is_weekend')
prophet.add_seasonality(name='daily_weekday',  
                        period=1,
                        fourier_order=5, 
                        condition_name='is_weekday')

# fitting the model
prophet.fit(train_prophet);

# part of the dataframe on which we want to make predictions
future = test_prophet.drop(['y'], axis=1)

# predicting values
forecast = prophet.predict(future)

# see https://github.com/facebook/prophet/issues/999 for the matplotlib_converts()
pd.plotting.register_matplotlib_converters()

# plotting the seasonality components found
prophet.plot_components(forecast)
Out[23]:
In [24]:
# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_prophet.ds, y=test_prophet.y,
                         mode='lines',
                         name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=forecast.ds, y=forecast.yhat,
                         mode='lines', 
                         name='Test - Prediction'))

# adjust layout
fig.update_traces(line=dict(width=0.5))
fig.update_layout(title='Prophet Forecast of Hourly Energy Demand',
                  xaxis_title='Date & Time (yyyy/mm/dd hh:MM)',
                  yaxis_title='Energy Demand [MW]')
fig.show()

# quantify accuracy
print(f'MAPE for Prophet\'s predictions: {mape(test_prophet.y, forecast.yhat)}')
MAPE for Prophet's predictions: 9.06%
In [25]:
# interval length
interval = 24 * 7

# intermediary variables for readability
x_true, y_true = test_prophet.iloc[:interval].ds, test_prophet.iloc[:interval].y
x_pred, y_pred = forecast.iloc[:interval].ds, forecast.iloc[:interval].yhat

# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_true, y=y_true,
                         mode='lines',
                         name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=x_pred, y=y_pred,
                         mode='lines', 
                         name='Test - Prediction'))

# adjust layout
fig.update_traces(line=dict(width=0.9))
fig.update_layout(title=f'Prophet Intra-Day Forecast of First {interval} Hours of Energy Demand',
                  xaxis_title='Date & Time (yyyy/mm/dd hh:MM)',
                  yaxis_title='Energy Demand [MW]')
fig.show()

# quantify accuracy
print(f'MAPE for interval of the first {interval} hours: {mape(y_true, y_pred)}')
MAPE for interval of the first 168 hours: 6.71%
In [26]:
# interval length
interval = -24 * 7

# intermediary variables for readability
x_true, y_true = test_prophet.iloc[interval:].ds, test_prophet.iloc[interval:].y
x_pred, y_pred = forecast.iloc[interval:].ds, forecast.iloc[interval:].yhat

# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_true, y=y_true,
                         mode='lines',
                         name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=x_pred, y=y_pred,
                         mode='lines', 
                         name='Test - Prediction'))

# adjust layout
fig.update_traces(line=dict(width=0.9))
fig.update_layout(title=f'Prophet Intra-Day Forecast of Last {abs(interval)} Hours of Energy Demand',
                  xaxis_title='Date & Time (yyyy/mm/dd hh:MM)',
                  yaxis_title='Energy Demand [MW]')
fig.show()

# quantify accuracy
print(f'MAPE for interval of the last {abs(interval)} hours: {mape(y_true, y_pred)}')
MAPE for interval of the last 168 hours: 11.89%

Prophet's MAPE: 9.06%

1st week forecast MAPE: 6.71%

Last wekk forecast MAPE: 11.89%

7 - Método III: Métodos de predicción simples. Implementación en R

Para la parte del trabajo realizada en R me he decantado por utilizar el paquete forecast de Rob Hyndman https://github.com/robjhyndman/forecast. La documentación está muy detallada, y además es posible encontrar numerosos ejemplos en su libro Forecasting: Principles and Practice.

Como ya vimos en la sección 3 de de las pruebas de raíz unitaria. El objeto base del paquete forecast es la clase ts incluida en el paquete básico de R. Forecast generaliza ts introduciendo la clase msts. El primero es un objeto time-series con frecuencia regular y estacionalidad única. Mientras que msts aunque también es un objeto de frecuencia regular (no permite fechas ni valores vacíos como un objeto xts), permite pasarle un vector con las frecuencias más importantes (diarias, semanales, anuales, etc). Esta flexibilidad veremos que puede marcar la diferencia a la hora de realizar predicciones en función de la serie temporal a tratar.

Para las siguientes secciones del trabajo he seguido un proceso incremental. Primero he analizado la serie con los métodos más sencillos (aunque bastante robustos y eficientes) y he ido aumentando la complejidad de los modelos utilizados (Regresión armónica dinámica con series de Fourier y finalmente TBATS).

En los siguientes puntos de la sección 7, mostraré los resultados obtenidos con estos métodos más sencillos.

7.1 - Mean method

Este método consiste simplemente en dar a todos los valores futuros, el valor medio de la serie.

>>> pjme_mean    <- meanf(train_ts, h=24*245)

# Plot
>>> autoplot(train_ts) +
  autolayer(test_ts, series="Test") +
  autolayer(pjme_mean, series="Mean", PI=FALSE) +
  ggtitle("Power Demand (MW) over time [2016 - 2018]") +
  xlab("Year") + ylab("Demand in MW") +
  guides(colour=guide_legend(title="Forecast"))

# MAPE (%) computation
>>> acc1 <- round(accuracy(pjme_mean, test_ts),2)[2,5]
>>> sprintf("Mean method MAPE: %s %%", acc1)

"Mean method MAPE: 14.7 %"

7.2 - Naïve method

En este caso, todas las predicciones reciben el valor del último dato disponible.

>>> pjme_naive   <- naive(train_ts, h=24*245)

# Plot
>>> autoplot(train_ts) +
  autolayer(test_ts, series="Test") +
  autolayer(pjme_naive, series="Naïve", PI=FALSE) +
  ggtitle("Power Demand (MW) over time [2016 - 2018]") +
  xlab("Year") + ylab("Demand in MW") +
  guides(colour=guide_legend(title="Forecast"))

# MAPE (%) computation
>>> acc2 <- round(accuracy(pjme_naive, test_ts),2)[2,5]
>>> sprintf("Naive method MAPE: %s %%", acc2)

"Naive method MAPE: 15.18 %"

7.3 - Naïve with drift method

En este método se permite aumentar o disminuir el valor de las predicciones en el tiempo. Es como si trazasemos una línea entre el primer punto histórico de la serie y el último y extrapolásemos con la misma pendiente a valores futuros.

>>> pjme_drift   <- rwf(train_ts, h=24*245, drift=TRUE)

# Plot
>>> autoplot(train_ts) +
  autolayer(test_ts, series="Test") +
  autolayer(pjme_drift, series="Drift", PI=FALSE) +
  ggtitle("Power Demand (MW) over time [2016 - 2018]") +
  xlab("Year") + ylab("Demand in MW") +
  guides(colour=guide_legend(title="Forecast"))

# MAPE (%) computation
>>> acc3 <- round(accuracy(pjme_drift, test_ts),2)[2,5]
>>> sprintf("Naive method with drift MAPE: %s %%", acc3)

"Naive method with drift MAPE: 15.04 %"

7.4 - Seasonal Naïve method

Los visualización de las gráficas dadas por este método parecen más complicadas de lo que en realidad son. Pero simplemente los puntos de las predicciones son los valores en estacionalidades anteriores. Es decir, si tenemos un valor y_t+1 = y_t en la estacionalidad pasada. Por ejemplo los puntos en Diciembre de 2019 tendrán el valor de puntos en Diciembre de 2018.

>>> pjme_s_naive <- snaive(train_ts, h=24*245)

# Plot
>>> autoplot(train_ts) +
  autolayer(test_ts, series="Test") +
  autolayer(pjme_s_naive, series="Seasonal naïve", PI=FALSE, alpha=0.7) +
  ggtitle("Power Demand (MW) over time [2016 - 2018]") +
  xlab("Year") + ylab("Demand in MW") +
  guides(colour=guide_legend(title="Forecast"))

# MAPE (%) computation
>>> acc4 <- round(accuracy(pjme_s_naive, test_ts),2)[2,5]
>>> sprintf("Seasonal Naive method MAPE: %s %%", acc4)

"Seasonal naïve method MAPE: 13.12 %"


8 - Método IV: Regresión armónica dinámica (DHR) para estacionalidades múltiples - Implementado en R

Cuando hay estacionalidades muy largas, una regresión armónica dinámica utilizando series de Fourier para modelizar los regresores externos es más indicado que métodos como ARIMA o ETS. En el dataset escogido para este trabajo, al tener datos por hora, estamos hablando de que cualquier estacionalidad presente va a ser muy alta, m=24 para estacionalidad diaria, m=168 para una estacionalidad semanal, m=8766 para estacionalidades anuales. La principal diferencia es que los métodos de la familia ARIMA o ETS están concebidos para estacionalidades más bajas (mensuales, m=12 o trimestrales m=4). La principal razón es de coste computacional, ya que estos métodos necesitan calcular m-1 parámetros, con el consiguiente coste computacional y de memoria.

En los métodos de regresión armónica, la componente estacional se calcula añadiendo transformadas de Fourier en las distintas frecuencias de cada estacionalidad modelada. El movimiento periódico de corta duración de la serie se modeliza con un modelo ARMA (modelo autoregresivo de media móvil).

Las principales ventajas que ofrece este método son:

  • Estacionalidad de cualquier duración (para diferentes estacionalidad, se añaden distintos términos de Fourier con distintas frecuencias).

  • La suavidad del patrón estacional puede controlarse con K, el número de pares de cosenos y senos en cada Transformada de Fourier para cada frecuencia. Un valor más bajo de K implica un modelado más suave.

  • El movimiento periódico a corto plazo se modela con un simple modelo ARMA.

8.1 - Implementación en R

Para implementar este método, me he decantado por utilizar la librería forecast() de Rob Hyndman. La forma de implementarlo es la siguiente:

Primero utilizamos la función auto.arima().

fit <- auto.arima(train_ts, seasonal=TRUE, lambda = 0, xreg = fourier(train_ts, K=c(12,2)))

xreg, es el parámetro más importante porque es el que se encarga de modelizar las estacionalidades y añadirlas al modelo ARMA.

La componente ARMA del modelo lo selecciona de forma automática utilizando como criterio la minimización del AICc (bias corrected Akaike information criterion).

>>> fit # Output resumen del modelo:
Series: train_ts 
Regression with ARIMA(5,1,5) errors 
Box Cox transformation: lambda= 0 

Coefficients:
         ar1      ar2     ar3      ar4     ar5      ma1      ma2      ma3    ma4     ma5
      1.6553  -1.2152  1.4254  -1.2693  0.3028  -0.6828  -0.0202  -0.8966  0.289  0.4407
s.e.  0.0142   0.0229  0.0197   0.0202  0.0126   0.0134   0.0142   0.0070  0.013  0.0096
      drift    S1-24    C1-24    S2-24    C2-24    S3-24   C3-24   S4-24    C4-24    S5-24
      0e+00  -0.1268  -0.0817  -0.0644  -0.0042  -0.0041  0.0026  0.0044  -0.0036  -0.0012
s.e.  1e-04   0.0042   0.0043   0.0007   0.0007   0.0004  0.0004  0.0003   0.0003   0.0002
        C5-24   S6-24   C6-24  S7-24  C7-24  S8-24  C8-24   S9-24  C9-24  S10-24  C10-24
      -0.0043  -1e-03  -7e-04  7e-04  7e-04      0  1e-04  -1e-04      0  -1e-04       0
s.e.   0.0002   1e-04   1e-04  1e-04  1e-04      0  0e+00   0e+00      0   0e+00       0
      S11-24  C11-24  C12-24   S1-168  C1-168  S2-168  C2-168  S1-8766  C1-8766  S2-8766
           0       0       0  -0.0439  0.0174  0.0211  0.0233   0.0191  -0.0330   0.1163
s.e.       0       0       0   0.0053  0.0053  0.0025  0.0026   0.3038   0.2747   0.1426
      C2-8766  S3-8766  C3-8766  S4-8766  C4-8766
       0.0976  -0.0064  -0.0053   0.0058   0.0081
s.e.   0.1387   0.0936   0.0935   0.0695   0.0706

sigma^2 estimated as 0.0002021:  log likelihood=47630.77
AIC=-95167.55   AICc=-95167.28   BIC=-94804.29

Una vez tenemos el modelo, el siguiente paso es hacer una gráfica de las predicciones frente a los valores de test, y comprobar que los residuos tengan una apariencia de ruido blanco:

>>> fit <- auto.arima(train_ts, seasonal=TRUE, lambda = 0, xreg = fourier(train_ts,                         K=c(12,2,4)))

# Plot:
fit %>% 
  forecast(xreg=fourier(train_ts, K=c(12,2,4), h=24*245), level = 99) %>% 
  autoplot(include=945*24, PI=FALSE) +
  autolayer(test_ts, series="Test", alpha=0.8) +
  ylab("Demand in MW") + xlab("Date")

checkresiduals(fit)

"Dynamic Regression with Fourier terms MAPE: 1.04 %"



A la hora de interpretar los residuos de este modelo, como con otras técnicas de regresión, lo principal es que los residuos tengan una distribución, normal, aleatoria o de ruido blanco como se conoce comúnmente en el campo de procesado de señales.

9 - Método V: TBATS

Alternativa a la Regresión armónica dinámica con términos de Fourier. Está desarrollada por De Livera, Hyndman, & Snyder (https://robjhyndman.com/publications/complex-seasonality/).

La función tbats() del paquete forecast de R está completamente automatizada, aunque el usuario tiene cierta libertad para seleccionar algunos parámetros de entrada. La elección de estos parámetros reduce el número de combinaciones posible y por lo tanto disminuye el número total de modelos a elegir, por lo que la búsqueda del mejor modelo será más eficiente.

Cabe destacar la posibilidad de pasar a la función el parámetro use.parallel=TRUE. Con esta opción seleccionada tbats() aprovechará la capacidad multi-hilo del procesador, con la consecuente reducción del tiempo de cálculo.

>>> train_ts %>% 
      tbats(use.trend = FALSE, use.arma.errors = TRUE,
            use.parallel = TRUE) -> fit2

>>> fit2
TBATS(1, {5,4}, -, {<24,1>, <168,1>, <8766,1>})

Call: tbats(y = ., use.trend = FALSE, use.arma.errors = TRUE, use.parallel = TRUE)

Parameters
  Alpha: 0.1098577
  Gamma-1 Values: 0.000275106 -4.054186e-05 -0.0003767173
  Gamma-2 Values: 0.0003092811 -0.0002484964 -0.0001945771
  AR coefficients: 2.419955 -2.72113 1.520097 -0.248297 -0.125298
  MA coefficients: -0.445064 0.253288 0.504738 0.15711


# Forecast and Plot:
>>> fc2 <- forecast(fit2, h=24*245, level = 80)
>>> autoplot(train_ts) +
      autolayer(test_ts, series="Test") +
      autolayer(fc2, series="TBATS", PI=FALSE, alpha=0.7) +
      ggtitle("Power Demand (MW) over time [2016 - 2018]") +
      xlab("Year") + ylab("Demand in MW") +
      guides(colour=guide_legend(title="Forecast"))

"TBATS MAPE: 1.22 %"



Vemos como el modelo TBATS ha tenido peor rendimiento que el modelo DHR. Como se puede ver en las imágenes no ha conseguido capturar de forma tan exacta la variación estacional de mayor amplitud.

9.1 - Resumen de los resultados

Model MAPE (%) Language
Holt-Winters 13.11 Python
Prophet 9.06 Python
Mean 14.7 R
Naïve 15.18 R
Naïve with drift 15.04 R
Seasonal naïve 13.12 R
Dynamic harmonic Reg. + Fourier 1.04 R
TBATS 1.22 R

10 - Conclusiones

Las conclusiones podría resumirlas en los siguientes puntos:

  • El modelo de Holt-Winters, es un modelo que no está pensado para tratar series temporales con estacionalidad múltiple, aún así ha parecido sacar una solución bastante robusta.

  • La librería Prophet de Facebook, ofrece de manera sencilla una doble API (en Python y R) y muy potente un modelo de predicción con variables de estacionalidad por Transformadas de Fourier. El modelo ha resultado ser muy robusto, sin apenas haber experimentado con los parámetros, y además ofrece una gran flexibilidad, añadir un gran número de estacionalidades de diferentes frecuencias, funciona con datos intra día automáticamente si se le inyecta la serie temporal en forma de dataframe pandas con un índice de tipo Datetime.

  • El paquete forecast de R otorga al usuario la posibilidad de trabajar de forma muy potente con gran cantidad de modelos, transformacioens y utilidades para analizar y predecir series temporales. Al igual que con la librería Prophet, la experiencia y el tiempo son factores clave para seleccionar los parámetros adecuados para cada serie temporal.

  • Se ha podido observar durante la realización de este trabajo la gran cantidad de factores y parámetros que hay que tener en cuenta para ajustar los modelos de predicción a una serie temporal. Por esta razón modelos como el de Facebook, cuyos valores han sido elegidos en base a ejemplos en la documentación, son con un gran alto de probabilidad subóptimos. Es muy notable la diferencia con los modelos automáticos implementados desde el paquete forecast de R. La selección automática a través de la minimización del AICc ha dado lugar a la obtención de modelos más indicados para la serie temporal que se estaba analizando.

11 - Futuras líneas de trabajo

  • En este proyecto se ha considerado una serie temporal univariante, sin ningún tipo de serie adicional variable complementaria ni covariantes (por ejemplo variables climatológicas como la temperatura, horas de sol, etc). Este es uno de los puntos débiles de este análisis, porque es precisamente la fuerte correlación entre el tiempo y el consumo eléctrico uno de los puntos fuertes de los modelos de demanda de electricidad. Este comportamiento se ha visto reflejado durante el EDA, cuando los patrones de gasto eléctrico diferían mucho entre las estaciones del año.

  • Los fines de semana, los periodos estivales (cuando es mayor la cantidad de gente de vacaciones), los eventos de un día afectan en gran medida al consumo. Así como la propia localización de los datos. No es lo mismo tomar datos de demanda eléctrica en una urbe, que en una zona rural, o en una zona turística donde se esperaran aún mayor varianza entre los periodos de vacaciones y el resto del año. Por todos estos factores es muy común la presencia de outliers presentes en los datos en forma de picos de gasto o días de muy poco consumo. Este es sin duda uno de los puntos fuerte de la Prohet de Facebook, que es la facilidad para modelizar este tipo de casos especiales.

  • Automatizar la selección del hiperparámetro K para los modelos en los que se modelizan regresores externos de estacionalidad mediante series de Fourier. Implementando un pequeño Grid Search entre diferentes valores de K, y tomando como condición una minimización del valor del AICc podría dar lugar a mejores predicciones.

  • Por último, pienso que sería interesante automatizar una validación cruzada con diferentes rangos de la serie temporal. De esta forma se podrían ejectuar diferentes modelos con diferentes rangos de datos y obtener finalmente un modelo más robusto.

La implementación de la validación cruzada podría partir de este ciclo junto con el fragmento de código que hay en el Anexo:

for train_index, test_index in tscv.split(X):
    print("TRAIN:", train_index, "TEST:", test_index)
    X_train, X_test = X.iloc[train_index], X.iloc[test_index] # pandas doesn't support directly integer based slicing, so that .iloc is mandatory
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

    X_train, y_train = model_forecast.fit()
    y_pred = model_forecast.predict(X_test) # in literature is usually seen y_pred as y_hat
    cross_validaiton_result(y_test, y_pred)

12 - Referencias

Anexo: Validación cruzada anidada para series temporales

En el análisis de series temporales es muy importante tener en cuenta la ordinalidad de los datos. Es precisamente ésta la propiedad característica de una serie temporal, la relación de cada punto con el todo el histórico de datos anterior a él.

Durante este proyecto, y por limitaciones de tiempo y capacidad de computación, he realizado las predicciones y los modelos con una división simple en un único conjunto de entrenamiento y un cojunto de validación. Sin embargo existe un procedimiento más sofisticado, que es la validación cruzada anidada, que podemos visualiar en la imagen siguiente:

En análisis de series temporales la proporción de datos incluída en los conjuntos de entrenamiento y validación suele ser del 80%-20%. Aunque suele haber desviaciones respecto estas cantidades dependiendo de la duración de la serie temporal, las características de estacionalidad y estacionareidad, variables covariantes, etc.

A continuación he dejado un pequeño fragmento de código que aprovecha la funcionalidad de la clase TimeSeriesSplit() del paquete scikit-learn para implementar una validación cruzada.

In [27]:
def plot_cv_indices(cv, X, y, ax, n_splits, lw=10, cmap_cv=plt.cm.coolwarm):
    """Create a sample plot for indices of a cross-validation object."""

    # Generate the training/testing visualizations for each CV split
    for ii, (tr, tt) in enumerate(cv.split(X=X)):
        # Fill in indices with the training/test groups
        indices = np.array([np.nan] * len(X))
        indices[tt] = 1
        indices[tr] = 0

        # Visualize the results
        ax.scatter(range(len(indices)), [ii + .5] * len(indices),
                   c=indices, marker='_', lw=lw, cmap=cmap_cv,
                   vmin=-.2, vmax=1.2)

    # Formatting
    yticklabels = list(range(1, n_splits+1))
    ax.set(yticks=np.arange(n_splits) + .5, yticklabels=yticklabels,
           xlabel='Sample index', ylabel="CV iteration",
           ylim=[n_splits, -.2])
    ax.set_title('{}'.format(type(cv).__name__), fontsize=15)
    ax.legend([Patch(color=cmap_cv(.8)), Patch(color=cmap_cv(.02))],
              ['Testing set', 'Training set'], loc=(1.02, .8))
    # Make the legend fit
    plt.tight_layout()
    
    return ax
In [28]:
# Let's separate the features_target dataframe in features matrix X and target vector y
def separate_features_target(df, target_labels):
    """
    This function takes the features+target DataFrame and splits in the features and target matrices
    """
    X = df.drop(target_labels, axis=1)
    y = df[target_labels]
    
    return X, y

X, y = separate_features_target(features_target, 'demand_in_MW')
display(X.head())
display(y.head())
date hour dayofweek dayofweek_string month quarter season season_string dayofyear dayofmonth weekofyear year
2002-01-01 01:00:00 2002-01-01 01:00:00 1 1 Tuesday 1 1 0 winter 1 1 1 2002
2002-01-01 02:00:00 2002-01-01 02:00:00 2 1 Tuesday 1 1 0 winter 1 1 1 2002
2002-01-01 03:00:00 2002-01-01 03:00:00 3 1 Tuesday 1 1 0 winter 1 1 1 2002
2002-01-01 04:00:00 2002-01-01 04:00:00 4 1 Tuesday 1 1 0 winter 1 1 1 2002
2002-01-01 05:00:00 2002-01-01 05:00:00 5 1 Tuesday 1 1 0 winter 1 1 1 2002
2002-01-01 01:00:00    30393.0
2002-01-01 02:00:00    29265.0
2002-01-01 03:00:00    28357.0
2002-01-01 04:00:00    27899.0
2002-01-01 05:00:00    28057.0
Freq: H, Name: demand_in_MW, dtype: float64
In [29]:
fig, ax = plt.subplots(figsize=(7,4))
n_splits = 5 #Number of train/cv/test folds
#max_train_size = int(len(features_target)*.67/n_splits)
tscv = TimeSeriesSplit(n_splits, max_train_size=None)

for train_index, test_index in tscv.split(X):
    print("TRAIN:", train_index, "TEST:", test_index)
    X_train, X_test = X.iloc[train_index], X.iloc[test_index] # pandas doesn't support directly integer based slicing, so that .iloc is mandatory
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]

plot_cv_indices(tscv, X, y, ax, n_splits)
TRAIN: [    0     1     2 ... 24229 24230 24231] TEST: [24232 24233 24234 ... 48461 48462 48463]
TRAIN: [    0     1     2 ... 48461 48462 48463] TEST: [48464 48465 48466 ... 72693 72694 72695]
TRAIN: [    0     1     2 ... 72693 72694 72695] TEST: [72696 72697 72698 ... 96925 96926 96927]
TRAIN: [    0     1     2 ... 96925 96926 96927] TEST: [ 96928  96929  96930 ... 121157 121158 121159]
TRAIN: [     0      1      2 ... 121157 121158 121159] TEST: [121160 121161 121162 ... 145389 145390 145391]
Out[29]:
<matplotlib.axes._subplots.AxesSubplot at 0x1fe838ff940>